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Magnetic flux trapped on the surface of superconducting rotors of the Gravity Probe B (GP-B) 
experiment produces some signal in the SQUID readout. For the needs of GP-B error analysis and 
simulation of data reduction, this signal is calculated and analyzed in the paper. We first solve a 
magnetostatic problem for a point source (fluxon) on the surface of a sphere, finding the closed form 
elementary expression for the corresponding Green's function. Second, we calculate the flux through 
"q^" ■ the pick-up loop as a function of the fluxon position. Next, the time dependence of a fluxon position, 

| caused by rotor motion according to a symmetric top model, and thus the time signature of the flux 

. are determined, and the spectrum of the trapped flux signal is analyzed. Finally, a multi-purpose 

program of trapped flux signal generation based on the above results is described, various examples 
of the signal obtained by means of this program are given, and their features are discussed. 



I. INTRODUCTION 

> ' 

The Gravity Probe B (GP-B) satellite is scheduled to fly in the year 2000. It contains a set of gyroscopes intended 
to test the predictions of general relativity that a gyroscope in a low (altitudes 650 km) circular polar orbit will 
precess, relative to a distant star, about 6.6 arcsec/year in the orbital plane (DeSitter, or geodetic, precession) and 
about 42 milliarcsec/year perpendicular to the orbital plane (Lense-Thirring, or frame-dragging, precession). To 
^\ ' provide the desired measurement accuracy (1 part in 10 5 for the geodetic effect), a magnetic London moment readout 
Q\ , using SQUID has been chosen, so that the experiment will be carried out at low temperature (~ 2.5° K), and the 
gyro rotors will be superconducting (see jlj], || for the design and status of the experiment; the history of GP-B 
development is found in [Q], and a survey of space relativity tests is in Q). The direction of the magnetic London 
moment developed in a rotating superconductor coincides with the direction of the rotation (spin) axis (F.London §; 
I for basic superconductor physics see 0; the description of gyromagnetic effects can be found in ||], Ch. 4). The 
■ corresponding magnetic flux through the pick-up loop of the SQIUD is proportional to the sine of the angle between 
the London moment vector and the pick-up loop plane, so the change of this angle, and thus the drift of the gyroscope 
• . axis, can be detected from the SQUID signal at the roll frequency of the spacecraft which will be deliberately rotated. 
However, along with the London moment dipole, there will also be quantum-size sources of magnetic field (fluxons) 
pinned to the surface of the superconducting rotor (see Ch. 5; ||, Ch. 12) which produce additional magnetic 
flux through the pick-up loop called trapped flux; its time signature will be present in the SQUID output. The low 
frequency part of this signal, though comparatively small under the GP-B conditions, might corrupt the accuracy of 
the London moment readout. On the other hand, its high frequency part can provide additional information significant 
for the experimental results. To make sure the trapped flux does not affect the measurement precision, as well as to 
extract useful information from it, one has to analyze the trapped flux signal and develop the code generating it, for 
the use in simulations of the GP-B error analysis and data reduction. This is the aim of the present paper. Note that 
the first work on the analysis of the trapped flux from a GP-B rotor was done by L.Wai in his thesis JTo[ ] . 

In sec. II we give a closed form solution to a magnetostatic problem of a fluxon on the surface of the gyroscope. In 
sec. Ill the solution is used to find the trapped flux signal in the pick-up loop as a function of the fluxon's position. 
The closed form expression for the trapped flux appears to be not very useful for further applications, so various 
exact and approximate formulas are also obtained. In sec. IV we investigate the motion of a fluxon with respect to 
the pick-up loop, thus finding the time signature of the trapped flux signal; we then go on to analyze its frequency 
spectrum. The last section contains a brief description of the program used to simulate trapped flux for the GP-B 
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data processing routines. Pictures of the high frequency signal, its low frequency envelope, and various Fourier spectra 
are presented and discussed. 



II. GREEN'S FUNCTION OF THE MAGNETOSTATIC PROBLEM 



The GP-B experiment will be conducted at low temperatures, so the fluxons can be treated as static (welded to 
the rotor's surface) and non-interacting ones. In such a case the total fluxon field is a superposition of the fields 
of individual fluxons. In addition, the rate of change of this field due to the rotor's motion is negligible, hence the 
magnetostatic approach should be used. Thus we consider a single fluxon whose characteristic size is on the order 
of 10 -5 cm ( ||, p. 184); due to a macroscopic size of the gyroscope (1.91 cm radius), the fluxon can be treated as 
a point source of magnetic field with the coordinate angles iff on the surface r = r g of the rotor. The spherical 
coordinates r, i9, ip here correspond to a Cartesian frame {x,y, z} fastened to the pick-up loop so that the origin 
coincides with the loop center and the z axis is perpendicular to the loop plane; the real relative motion of the fluxon 
and the loop, i. e., the dependence of the fluxon position angles ipf on time, will be incorporated and examined 
in sec. IV. 

In these settings, the boundary value problem for the magnetic potential W(r) of the fluxon outside the rotor is 
formulated as 



A* (r) = 0, 



dr 



T > Tg, < Of < it, < iff < 2tt 
$ 



r^sim? f 



5(0-0f)8(<p-<pf), 



where $o = h/2e is the magnetic flux quantum, and the magnetic field is 

B = -V# 



(1) 
(2) 

(3) 



Evidently, up to a factor $o, ^ is the Green's function of the external Neumann boundary value problem for a sphere. 
A standard separation of variables leads to the following series representation of the solution to (flh, (0): 



(r) = tf (r, 0, ip) = —5- E ( Mlm cos m ^ + Nlm sin m ^ (—) P ™ ( cos § ) 



1=0 m=0 



with the coefficients given by 

21 + 1 (l-m)\ 



(1 + tfmo) {I + 1) {l + m)\ 



(costf / ) cos rrupf , N tr 



21 + 1 (l-m)\ 
(l + l) {l + m)\ 



Pn (cos$f)smm(p f 



(4) 



(5) 



As it turns out, this series may be summed to give the closed form solution for ty. To determine it, we first introduce 
3) into (0) to obtain 



* (r) 



■E 



21 + 1 
l + l 



Pi (cos i?) Pi (cos $ f ) + 2j2 p r ( c ° s p r ( c ° s fif ) c ° s m ( ^ - w ) 



m=0 



Then, by applying the addition theorem for Legendre functions (see O], 10.11, (47)), we convert the latter into 



* (r) = > 



1 



1=0 



l + l 



l+l 



JL ) P;(cos 7 ) 



_$o_ 

47rr„ 



L (=0 1=0 



where 7 is the angle between the directions to the fluxon and to the observer: 

COS7 = costfcosi?/ + sim?sim9/cos (if — iff) 



(6) 



(7) 



The first of the series in the above expression for 'J is obviously the generating function for Legendre polynomials 
(see [0, 10.10, (39)), the second one is just an integral of it, namely, 



In 



i-c 
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Using these results in 
elementary functions: 



we can now write the magnetic potential in its final form as a finite combination of 



(r) = $ G(r,r f ) 



2tt 



r - r f 



— ln-5- 

2r 



r ■ r f + r g |r - r f | 
rr„ — r • rf 



(8) 



where G(r, rf) is the mentioned Green function and rf = {r g ,i}f, iff} is the position vector of the source. The first 
term here, as one would expect, is a half of the potential of a point charge, and the addition to it describes the 
contribution of the curved boundary. 

Since, surprisingly enough, we were not able to find this explicit formula in literature, it seems reasonable to give 
here a closed form expression for the Green function of the corresponding Dirichlet problem (Gd), in which the 
boundary condition (0) is replaced by 



4>o 



r g sin ■d 



-5{0-0f)8(<p-<pf) 



The result then is 



*(r) = $ G D (r,r { ) 



47T 



(9) 



(10) 



Note that Green's functions for the corresponding internal problems can be obtained from (JSJ) and (JlOj) by means of 
inversion. 



III. TRAPPED FLUX AS A FUNCTION OF A FLUXON POSITION 



Magnetic flux measured by the pick-up loop of a GP-B SQUID is the flux through the circle of the radius R in 
the plane z = 0, or, equivalently, the flux through the (upper) hemisphere. The dependence of the trapped flux on 
the fluxon position turns out to be rather complicated, especially for the GP-B design, when the gap between the 
rotor and the loop is very small as compared to the pick-up loop radius R. For that reason we give here a number of 
different representations of the trapped flux as a function of the fluxon position; each of them has its own merits and 
drawbacks and is thus used for different purposes pertinent to our investigation. 



A. Trapped flux in terms of series of Legendre polynomials 



The simplest way to calculate the trapped flux is to integrate over the hemisphere the series expression for the 
radial component of the magnetic field obtained from (H)-(||): 



B, 



dA = 



hemisphere (r—R) 



r=R 



9* 
dr 



hemisphere (r—R) 



r=R 



oo i /"l 

dA=$ ^(l+l)(&\ M m / Pi( a )da; 
i=o ■'° 



all spherical harmonics with m ^ here have averaged out over the azimuthal angle if. The last integral is calculated 
with the help of the known relations of the theory of Legendre polynomials (see pT| , 10.10, (14), (2), (4)) : 

p ,„x Pl+iW - PLi(s) pm 1 p tn , n p rn s (-i) fc r(fc + i/2) _ 

Pi{s) = r— ; Pi(i-) = 1 P 2 fc+i(0) = 0; P 2 fc(0) = — -= n ; <,fc = 0,l,... 

I + 1 y/ir k\ 

r(C) is the Euler gamma-function. Then, after inserting the values Miq from (|^), we arrive at the following expressions: 

<f>„ 



$ / (cosi? / ) = ^F <5 (cos l ? / ); 



F 5 (s) = £ (1 - <5) 2fe+1 P 2k+1 (s) [P 2k (0) P 2fc+2 (0)] = -= ]T(-1) 



k ^^nk+l/2)(l-6f k+1 P 2k+1 (s) (11) 



k=0 v k=0 

Here S denotes the dimensionless gap between the pick-up loop and the rotor, < 6 = (R — r g )/R < 1. 
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From the point of view of signal processing, Fg(s) is a transfer function which converts the "input" fluxon position 
signal Si n (t) = cos d / (t) (the position is changing with the time as the rotor moves relative to the pick-up loop, see the 
next section), into an "output" trapped flux signal S out (t) — 0.5&oFs(Si n (t)) which is present in the GP-B readout. 
For the reason that the total contribution to the flux of any number of fluxons scattered in any way over the rotor's 
surface is given by the sum of the values of the same function Fg taken at proper different values of its argument, 
it was called "universal curve" in [pLOf] - Clearly, Fg(s) is an odd function of s; in particular, Fg(0) = means that a 
fluxon sitting exactly in the pick-up loop plane does not, of course, register any flux. 

By setting 6 = in (|ll]) (the loop on the surface of the rotor), we immediately find 

2 °° A- + 3/4 C 1 if < s < 1; 

W^EHf^ + l/^). if S = 0; (12) 
v fc=o v ' \ — 1 if — 1 < s < 0. 

(the last equality here is proved by expanding its right-hand side in orthogonal series of Legendre polynomials). 

This result obtained by L.Wai JTu] has a clear physical meaning: when the pick-up loop lies on the rotor's surface, 
same as the point source of field always does, the flux through the loop remains unchanged (±$ /2, half of the total) 
while the fluxon stays in either of the hemispheres separated by the plane of the loop, and changes it sign by a jump 
when the fluxon crosses this plane. However, equation ( fl2"| ) also demonstrates the difficulties in using expression jll[ ) 
for GP-B, where 6 = 0.025 is very small: for any S > the series ( |Tl"| ) has an absolutely converging majorant, 
so its sum Fg(s) is an analytical function of s, but it has a jump discontinuity at s = when (5 = 0. Therefore 
the series ([ll]) converges worse and worse with the separation 6 becoming smaller and smaller, which makes ( |ll| ) 
practically unacceptable for accurate numerical calculations at the required value of separation. It also turns finding 
a uniform in s asymptotic expansion of F$(s) for S — > into a rather difficult mathematical problem. The effect is 
that for small positive values of 6 the transfer function has a shape of a very steep "kink" (recall that Fg (s) is odd): it 
is almost constant outside a small vicinity (—As, As) of the origin, with A,5 = 0(6) as shown below, and is equal to 
zero at s = with a huge gradient ~ 0(1/6) there (see fig. 1). That is why we are deriving three more representations 
for Fs(s) in the following subsections. 



B. Integral representation of the trapped flux 



An integral expression for Fg(s) is obtained by replacing the Legendre polynomials in ( p"l| ) by their integral repre- 
sentation (see [In), 10.10, (43)) 



-P2fc+i(cos $/) 



1 exp[i(2fc+ 1 + 1/2)] dip 

7T J-tff y/2 (cos?/; - cos $ f) 



Changing then the order of summation and integration, we arrive at a sum of two hypcrgcomctric scries which are 
readily summed up to result in: 



F s (cos# f ) 



$oV2 



dtp exp (iip/2) 
o \J cos tp — cos t9 / 



A 



Vi + P 



VTTA2 

2A 



1 

2A 



A = (1 — 6) exp (iip) 



(13) 



Representation (O) is very convenient for precise numerical calculation (and, in fact, is used for this purpose in 
our code, see sec. V), because the integrand in ( |l3| ) is an algebraic one, and the weak singularity at the upper limit 
can be taken care of rather easily. The plot of the transfer function computed from (|l3|) is given in fig. 1, along 
with the graphs of its various approximations described in the next subsection. The relative error of the numerical 
computation has been kept within 10 -5 . 



C. Elementary approximations of the trapped flux 

From the described behavior of Fs(s) for small 5 it is clear that to effectively approximate it one needs the value of 
its gradient at s = and the "saturation" value ^(1), in the first place. Fortunately, it is possible to compute these 
quantities exactly, and they are 



fs = F s (l) = 



1 



1 - 



26 -6 2 



= 1 - (V2- l)6 + 0(6 2 



(14) 



4 



KS 



dF 5 (s) 



l^| E(1 -,-K (1 -, 



2 
n 



2 + 0(SlogS- 1 ) 



0; 



(15) 



here K(fc), E(fc) are complete elliptic integrals of the first and second kind, respectively (see JL4|, Ch. IX for their 
definitions and asymptotic behavior at k — > 1 — 0). The formulas are derived from ([nj) by the direct summation of 
the corresponding series of Legendre polynomials carried out in the Appendix. 

The simplest approximation of the transfer function for 5 — > +0 is evidently a piecewise-linear one, 




with As defined in a natural way as 



ks A s = fs, 



if A s < s < 1; 

if \s\ < A s ; 

if -1 < s < -As, 



A s = ^ = ls + 0(d 2 ) 
ks 2 



(16) 



(17) 



It turns out that this approximation gives the right qualitative picture of the signal and even is not too bad quantita- 
tively, providing, for all values \s\ < 1, the error within 1/3 for both S — 0.3 and 5 — 0.025. This accuracy, however, 
is not enough for the GP-B simulations, moreover, the largest error, associated with the jump of the derivative of 
function (|16| ) at s = ±As, occurs in a very sensitive transition region where the fast growth of Fs(s) is replaced by 
its almost constant behavior. 

A much more attractive approximation is given by the function 



zr / \ f , (kkss 
Fs[s) w -f s arctan - — 

7T V 2 IS 



-0 



(18) 



The parameters here are arranged in such a way that the slope at s — is exactly ks and, in the spirit of asymptotic 
methods, the true saturation value is achieved when ks-s = oo (note that another "simple and natural" approximating 
function, the hyperbolic tangent, is not acceptable, because the rate of approaching of fs by Fs(s) is a power rather 
than exponential one). The performance of the approximation (18) exceeds all expectations, giving, over the whole 
range of s, the maximum error of 20% for S = 0.3, and only 1.8% for 6 = 0.025. 
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Universal Curve Fs(s). 



The accuracy is mostly lost outside the transition zone (—As, As) due to the fact that fs is achieved only at infinity. 
This can be dealt with by redefining the parameters to have both the exact slope at s = and the right value at 
s = 1, which produces 



5 



Fg(s) ~ ^4,5 arctan 



ngs 
~As~' 



Ag arctan 



us 



fs 



-0 



(19) 



This 'adjusted' arctan gives the maximum error within 0.3% for 5 = 0.025; and even for as large a separation as 
5 = 0.3 the error is still about 0.6%. Same as ( |l6| ) and (|l8|), the dependence ( |l9| ) is shown in fig. 1. As versus 6 is 
plotted in fig. 2; note a relative flatness of the of the function. 




Fig. 2. Dependence of As on 6. 



D. Closed form expression of the trapped flux 



The explicit formula for the trapped flux can also be obtained, though not that easily, from equation (|l l|) , however, 
a direct way to get it is to integrate the closed form expression for the magnetic field through the pick-up loop plane 
z = 0. For this plane $ = 7r/2, r — p (the polar radius); in addition, we can redefine tp by setting tp* — 0. Then 
equations (JsJ) and (^J) provide the needed component of the magnetic field in the form: 



B, 



z=0 



$0^s cos$/ 
2^ 



1 



, sin § f cos tp 



sin i? f cos tp 



1 



X 3 (p,tp) 2rlpX(p,v)Y + (<p)Y-(<p) 2r 2 q pY+(<p)Y-(<p) 2r 2 qP Y-(<p) 



where 



X(p, tp) — ^jr g — 2r g psm'df cos(p + p 2 , Y± (tp) — 1 ± sin §f cos cp 



(20) 



(21) 



Now we need to integrate (f20|) over the area of the pick-up loop. First we calculate the simple, though rather 
cumbersome, algebraic integral of the field ( pp[ ) times pdp over the polar radius from to R (if instead one first 
integrates over tp, elliptic integrals of a complicated argument appear in the result that make the closed form radial 
integration very difficult). As we are then to integrate over the period of cos tp, the terms odd in cos tp can be omitted, 
and we obtain: 



$/(cos 1 ? / ) = ^F 5 (cos^ / ) 



§0 r g C0S1?/ 

2^ 



2 IT 



dtp 



R 2 



R 



2r 2 X(R,tp)Y + ( V )Y-(tp) 2r 2 Y.(tp) 



In view of (|T]), this integration is also rather straightforward and leads to the desired result: 

Socostf/f 1 2S-S 2 \U(v + ,k) : II(i/_,fc) 



*/(cosi?/) 



2 1-6 l|cos0/| 7 r x /2(l- ( 5)(l+sin?9 / ) + l 5 2 



1 + sin <&f 1 — sin ■& t 



(22) 



(23) 



G 



where 



. 2sini?/ / 4(1 - 8) sin??/ 

^ ( ^ ) = T T±^4' k ^ = p(i-S)(i + ^) + P (24) 

and II(^, fc) is the complete elliptic integral of the third kind (see |l4[ , Ch.IX). As a consistency check, one may 
calculate the saturation value and the derivative at zero of the transfer function rt23) to see that they are indeed equal 



to the previously obtained values (14) and ( |15| ) 

The first term in fl23|) evidently has a jump discontinuity at s = cosi?/ = 0. Therefore, for all finite S, the second 
term must contain the discontinuity of the opposite sign, to make the sum of two analytical in s. Hence for small 
positive S in the transition zone we are dealing with a small difference of two large quantities, which is always a 
problem. Also, the first term in ( f23|) coincides exactly with the expression (12) for 6 = 0, hence the second one should 
disappear in this limit, which it necessarily does in a very nonuniform way. Evidently, such an expression cannot be 
effectively used for both numerical and analytical purposes when 6 is small enough, which is our case. 



IV. FLUXON KINEMATICS AND SPECTRAL DECOMPOSITION OF THE TRAPPED FLUX SIGNAL 

Now we need to determine the time signature i? / (i) of a fluxon polar angle in the pick-up loop frame, to complete 
the investigation of the trapped flux signal. 

In doing that we use four Cartesian coordinate systems. The first one, {x, y, z}, has been introduced in sec. I; it 
is fastened to the pick-up loop, and z is the unit vector normal to the loop plane. The second coordinate system, 
{x r , y r , z r }, is associated with the roll axis of the spacecraft, Cj r — z r (fig. 3). The roll axis is almost in the pick-up 
loop plane, that is, the roll axis — pick-up loop plane misalignment a < 10 -5 is very small. The third set of coordinates, 
{xl , Vl, zl}, is related to the angular momentum vector L in a way that zl = L/|L|. Both r- and L-coordinates 
are fixed in the inertial space, since the roll axis is pointed to the Guide Star, and we can so far neglect the pointing 
errors, as well as the relativistic drift of L. We choose axes y r and yx in the plane containing both z r and Zl, then 
the perpendicular to this plane axes x r and coincide (fig. 4), and the following relations are true: 

zl ■ z r = y L ■ y r = cos/3 , z L ■ y r = -y L ■ z r = sin/? 

x r • z r = x r • z L = x r • yy = x r • y L = (25) 

Here (3q is the roll axis — angular momentum misalignment which is required to be < 5 x 10~ 5 rad in the GP-B 
experiment. 





Fig. 3. Mutual Orientation of Roll and Loop Coordinates. 
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Fig. 4. Mutual Orientation of Roll and Angular Momentum Coordinates. 

A symmetric top with the moment of inertia / + AI relative to the body symmetry axis and equal and slightly 
different value / for the moments of inertia about the other two axes is a very good model for the GP-B rotors (note 
that \AI\/I < 10~ 5 for them). Therefore, we choose the fourth Cartesian coordinate system {xb, Ub, zb} fixed in 
the rotor's body with zb directed along the rotor's symmetry axis. 

The dynamics of a symmetric rotor is well known and relatively simple (c. f. [|l2|]l3| l). Its motion in the L-coordinatcs 
is a precession about zl with the spin frequency 

u s = j, (26) 

and rotation about the rotor symmetry axis zb with the frequency 

L / AI 



Urot 



I + AI 



cos 7s ~ oj s (l — ) cos 7s; (27) 



< 7b < 7r is the angle between zl and zb- 

For the signal of the trapped field we need, however, the time dependence of the position of a fluxon in the inertial 
coordinates, hence we need expressions of xs(£), ys(i), zs(t) in terms of X£, y_L, ^l- The latter is found with the 
help of the Euler angles (see for instance [^2[ ) in the form 

z B {t) = z L cos 7b + x L sin 7s cos0 s + yl sin 7b sin0 s 

ys(*) = z l sin7BCOs0 p + 

Xi (cos 7b cos0 s cosOp — sin0 s sin0 p ) + y^ (cos 7b sm.9 s cos6 p + cosO s sin# p ) 

Xs(t) = — zl sin 7b sm9 p + 
Xi (cos 7b cos 9 S sin p + sin 9 S cos P ) + yL (cos 7b sin 9 S sin p — cos 9 S cos p ) (28) 
Here the spin and polhode phases are 

9 s (t) = u s t + 0°, p (i) =w y i + 0°, 0% tp = const, (29) 
and w p is a polhode frequency, 
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L \AI\ \AI\ 

UJ p = j —j— cos IB = w s — — COS 7 B 



(30) 



(In the body-fixed frame the instant angular velocity vector rotates around the rotor's symmetry axis with the polhodc 
frequency). Using this, we obtain the following expression for the unit vector e/ in the direction of a fluxon (i. e., 
of an arbitrary fixed point of the rotor surface at some polar, < £ < tt, and azimuthal, < r\ < 2tt, angles in the 
body-fixed spherical coordinates): 

et = z B (t) cos£ + (x s (i) cost? + y B (t) sin 77) sin£ = ei(t)x L + e 2 (t)y L + e 3 (t)z L 



e\{t) = sin£ cos7s cos9 s {t) su\(9 p (t) + ?/) + sin 9 s (t) cos(9 p (t) + 77) + cos£sin7B cos8 s (t) 



e2(t) = sin£ cos7b sm9 s (t) sm(0 p (t) + 77) + cos8 s (t) cos(8 p (t) + 77) + cos£sin7B sm9 s (t) 
e 3(t) = ~ sin^sin7B sin(6* p (t) + 77) + cos£cos7b 



(31) 



According to the results of sec. Ill, we only need the cosine of the angle $/(i) between ef(t) and the normal z(t) 
to the pick-up loop plane to study the trapped field signal; together with the loop, z(t) rotates about co r with the 
frequency lu t (see fig. S2): 

z(i) = cos(7r/2 — a)uj r + sin(7r/2 — a)(cos# r x r + sin(9 r y r ) = sinaz r + cosa (cos0 r x r + sin^y,.) 

9 r = 9 r (t) = u r t = roll phase (32) 

By means of this, ( |3l| ) and formulas (^5|) relating the r- and L-coordinates, to the first order in the misalignments 
/3o and a we obtain (quadratic and higher order terms are several orders below the required GP-B accuracy): 



cos7?/(t) = a s _ r sin0 s _ r + a(/3o sin# r + a) 

Q s - r (t) = (iU s — LO r )t + q s ~r] 9 r (t) — LU r t 



(33) 



For a perfectly spherical rotor AI = and the amplitudes and initial phase here are true constants whose values 
depend only on the position of a fluxon relative to the symmetry axis, a s _ r = sin£, q s - r = T], cl = cos£. If, on 
the other hand, AI 7^ 0, they start to vary slowly with the time at the polhode frequency according to 

a s - r (uj p t) = \J [cos £ sin j B + sin £ cos(w p i + 9° + 77)] 2 + sin 2 £ cos 2 73 sin 2 (ui p t + p +rj) 



tan q s _ r (uj p t) 



si n 



£ cos 7b sm(uj p t + 9° + ?;) 



cos£ sin7B + sin£cos(u;pi + 9 p + 77) 



a{u) p t) = cos £ cos 7b — sin £ sin 7b sin(o;pt + 9® + 77) 



(34) 



Note that under the conditions of the GP-B experiment the spin frequency is always much larger than the roll and 
polhode ones, cu r <~ 5 x 10~ 5 w s , uj p ~ 10~ 5 w s . Since generally the second term in the first of equations (21) is about five 
orders of magnitude smaller than the first one, the input signal for the trapped flux output $/ (i) = (<&o/2) Fg (cos $f(t)) 
is a single carrier harmonics of the (high) spin minus roll frequency (O s „ r ), slowly modulated in the phase and 
amplitude at polhode frequency, added by a small D.C. offset (era), and a small low frequency harmonics (9 r ), both 
modulated at lu p . Therefore it is natural and convenient to represent <&/(£) as a Fourier series of spin minus roll 
harmonics with the amplitudes modulated by low frequencies, namely: 

*/(*) = ^-P*(«* */(*)) = 



£0 
2 



(oj p t) 2_ \ Ak(uj p t) sin(2fc + l)Q s - r (t) + a{oj p t) (/3o sinw r i + a) 2_ \ Bk(u> p t) cos2kQ s ~ r (t) 



fc=0 



fc=0 
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A k (u p t) = 7T ^ 2k + 1 ^ I cos{2k + l)ip cosip Fg(a s _ r (uj p t) simp) dip + 0(/3%); 



B k (uj p t) = ——— cos2ki>F^a s - r (uj p t)smiP)diP + O([3 l ); (35) 
tt(1 + Oko) Jo 

here prime denotes the derivative of Fg(s) in s. 

As readily seen, the amplitudes of odd harmonics of s _ r (A k ) are generally of the order of unity and decrease as 
0(k~ 2 ) for the large enough number k. In contrast with that, the amplitudes of even harmonics, which are linear in the 
misalignments, are at least four orders of magnitude smaller but decrease only as 0(k~ 1 ), k — » oo. In addition, the even 
harmonics are modulated also by the roll frequency u ri so that, along with the harmonics 2fc0 s _ r (i), k = 0, 1, ... , 
with amplitudes a a r (uj p t) Bk(uj p t), harmonics 2fc0 s _ r (i) ± u) r t are present, whose amplitudes differ only by the 
misalignment involved, 0.5/?o instead of a. 

With all this in mind, one can easily understand that the full spectrum of the trapped flux signal consists of the 

following series of frequencies: (2fc+l)(o; s — w r )imw p , 2i(u s -u r )±w r ±mw p and 2fc(w s -u r )±mw p , m, k = 0, 1, 

The highest peaks are at (2k + l)(u) s — w r ), and those at 2k(uj s — uj r ) ± uj r and 2k(uj s — uj r ) are four to five orders of 
magnitude smaller. All of them are surrounded by an appropriately scaled forest of side bands separated by imw p . 

The only remaining thing is to discuss briefly the total flux $ produced by all fluxons. There are always some 
N pairs of fluxons and antifluxons present on the rotor's surface after cooling the rotor down below the transition 
temperature (the antifluxon is a fluxon with the opposite sign of the magnetic field) . Experiments have indicated that 
the expected number of the pairs is around N ~ 100, at the most. We denote any values related to either fluxons or 
antifluxons by indices / and a, respectively, numbering them with the index i — 1, 2, . . . , N; for instance, their body 
coordinates will be £j, rff and £*, rf a , the input signals S}(t) = cos?9j(t), S l a (t) — cos^(t), etc. 

The general expression for the total trapped field flux is given by 

N N 

*(*)=£[*/(*) + **-(*)] = ^^[F <5 (cos4(t))-F <5 (cos<(t))]; (36) 

i=0 i=0 

obviously, the full spectral representation of $(t) is just a scaled up version of $/(£) given in (jslf). 

Since for small 8 the transfer function Fg(s) is close to ±Fg(l) ~ ±1 everywhere except a small vicinity of the origin 
(see sec. Ill), expressions (|3^), ( [35] ) demonstrate that the maximum value of <I>(t) is distributed according to the usual 
counting statistics, provided that the distribution of fluxons over the surface of the rotor is the uniform random one. 
Therefore N fluxon-antifluxon pairs in this case should produce the total flux on the order of y/N$>a for 'large' N. 



V. CODE AND SIGNAL ANALYSIS 



For the GP-B error analysis and data reduction one needs to simulate the trapped flux signal as expected in the 
SQUID output. To do that, the results obtained in the previous sections were utilized for writing a program able to 
fast enough generate, store, and analyze the high-frequency signal. The code written in the MatLab v. 5.0, to ensure 
compatibility with other GP-B software, is available from the authors. 

The program is very versatile, allowing for many options and many different tasks. For instance, there may be a 
different number of fluxons, their positions may be read either from a prewritten file or generated at random according 
to different probability distributions. Transfer function may be calculated by means of several different expressions 
introduced in sec. III. Generation of the high frequency signal and/or its slow varying Fourier amplitudes (|35|), ( |3"6| ) 
is possible. In addition, all gyroscope and pick-up loop parameters (radii, rotor asphericity, misalignments, etc.), as 
well as the discretization frequency, time intervals, and all angular velocities may be specified in an arbitrary way. 

A lot of attention in the program's realization has been paid to the fact that tracing positions of as much as 200 
fluxons for long enough periods of time with high discretization frequency easily becomes too memory consuming. 
The program thus has been optimized in several directions, such as not to cause excessive memory swaps to the hard 
drive, not to lead to the memory fragmentation, and to access the hard drive for data storage as infrequently as 
possible. The following data may be useful to estimate the code's speed: on a Sun UltraSparc 5 with 128 Megabytes 
of RAM running System V, Rel. 4-0 and having a network mounted storage drive it takes, depending on the network 
load, 1.5 up to 2 hours to generate one hour of signal of 100 fluxon pairs at a sampling frequency of 2200-ffz (the 
actual sampling rate of GP-B electronics). 

Here we will not elaborate more on the code details but continue with the results of our simulations. All of them 
have been performed with the parameters set at the values expected for the GP-B experiment (see c. f. (!]-§)■ ln 
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particular, the spin frequency f sp m = 100 Hz, the roll period T r = 3 min, the polhode period T p « 43.6 win; recall 
that S = 0.025. 

In fig. 5 the signals are seen generated by different number of fluxons distributed in various ways over the surface of 
the gyroscope. In all of the graphs the 'adjusted arctangent' approximation (|l9| ) to the universal curve is used. Fig. 
5, a shows signals of a single fiuxon (without an antifluxon counterpart) positioned at different points on the gyro. 
The majority of fiuxon positions provide signals like the one drawn in the solid line in the figure. The dashed and 
dash-dotted lines correspond to rare fluxons oscillating in the small (~ As) vicinity of the pick-up loop plane, which 
is why their amplitude is smaller. On the average, one cannot expect too many fluxons like that, however, each of the 
four GP-B rotors will carry just one particular realization of the fiuxon position distribution, so these 'weak' fluxons 
are possible. 

Fig. 5,b shows various signals from one fluxon-antifluxon pair. Again, the solid line correspond to 'the most 
probable' signal: fiuxon and antifluxon are far from each other (though not opposite on the sphere) and have large 
oscillation amplitudes. 

Fig. 5,c shows typical signals of 5, 15, and 100 pairs distributed randomly with the uniform probability over the 
gyro surface. The VN growth of the signal is visible; the complexity of the signal profile also clearly increases with 
N. 



a. Single Fiuxon Signals 
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c. Many-Pairs Signals 
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d. GP-B Signal at Different Times 
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Fig. 5. Simulated Readout Signals. 
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Fig. 5,d shows short fragments of the 12 hours of signal generated for the test of the GP-B data reduction algorithms. 
There are 100 fiuxon pairs distributed unevenly: 60 of them are uniformly spread at random over the surface (just like 
in fig. 5,c), while the remaining 40 are used to create a total net flux of ~ 40 $o along some random axis. This should 
account for a small residual magnetization of the rotor at the time when it was made superconducting (see ]l7t). This 
magnetization not only significantly increases the amplitude of the signal, but also smoothes it out. Different curves 
in the figure correspond to the signals taken at different stages of the polhoidal motion (namely, 0, 15, and 24 minutes 
from some reference point) for a duration of 3 spin periods. 

In fig. 6 a low-frequency envelope is plotted of the signal from fig. 5,d used in GP-B simulations. The graph was 
constructed by splitting the magnetic flux signal into two-second blocks (4400 data points in each) and plotting the 
maximum value of the flux for each block. Periodicity of the large scale structures of the envelope with approximately 
the polhode period of about 43 min is clear. On the other hand, a comparison of the signal in any two corresponding 
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regions demonstrates that the short scale features, presumably introduced by the roll frequency and other less intensive 
harmonics, are not repeated precisely every polhode period T pi which is expected because T p and the roll period T r 
are incommensurable. 
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Fig. 6. Envelope of The Simulated Trapped Flux Signal, T p « 43.6 mm. 
a. Odd Harmonics b. Even Harmonics 
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Fig. 7. Slowly Varying Amplitudes of Fourier Harmonics of Trapped Flux Signal, T p « 43.6 min. 

Fig. 7 shows the slow polhoidal variation of Fourier amplitudes of the spin minus roll harmonics calculated according 
to (pq) and summed over the fluxons and antifluxons. The first ten odd and even harmonics are shown in plots a and 
b, respectively. Recall that in the expression (35) for the flux all even harmonics are multiplied by the misalignments, 
so that the actual vertical scale in fig. 7,b is about 10 5 of that in fig. 7, a. The pictures clearly show that the odd 
harmonics drop much faster with the number than the even ones, as predicted. It is interesting to note that the lowest 
even (n = 0) harmonics, which gives the amplitude of the D.C. and the roll frequency components, has a shape rather 
distinctive from the profile of the other modes. 
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APPENDIX. SUMMATION OF CERTAIN SERIES OF LEGENDRE POLYNOMIALS 



Here we give a derivation of formulas (|14|), (|15|) for f$ — Fg(l) and for the slope Kg of the transfer function at s = 0. 
We use the Pochgammer symbol (a)o = 1, ( a )k = a ( a + 1) . • • (a + As — 1) = T(a + k)/T(a), as well as the standard 
notation 



F(a, b, c;()=J2 



(o) fc (6fc) C k 



for the Gauss hypergeometric function of the argument £ and parameters a, b, c. From ( [ll|) we have 

F 5 ( S ) = F«( S )-if>( S ); 



fe=0 '"' v_/ fe k=0 y r \ / k 

where we introduced r) = 1 — S for brevity. 
Calculation of fg . Since -P«(l) = 1, we have 



Fi 1) (l) = 2r ? ^. 



and for the elementary expression of the hyper geo metric function we have used formula (11) from |D| , 2.11. with 
a = 1/2, 6=1. Combining these results with (|37|), we obtain 



2 



A - «(.) - F» W - Ff >(,) . -JL - £±f^l . i (i - -^L= 

which, in view of 77 = 1 — is exactly the expression (p"4|). 
Calculation of n s ■ As (see |0|, 10.10, (12)) 

(-l) fc /3 



P^ +1 (0) = (2fc + l)P 2fc (0) 



jfc! V 2 



A' 



from (|37h we find: 



dF s 



ds 



=0 fc=0 K ' w ^ fc=o K - 



,J, _ 2 V T?fi /o 1/01. „2n _ 4? 7 



2^(1/2, 3/2, 1; rf) = F(l/2, -1/2, 1; rf) = 1 Efo), (38) 
\ — r\ l 7r(l — 77^) 

where E(?j) is the complete elliptic integral of the second kind, and we have exploited the classical relation (see [ jl5| , 
2.1.4, (23)) 

f(o, 6, c; C) = (l-C) C ~ a ~^(c-a, c-6, c; C), 
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and the expression for the elliptical integral in terms of the hypergeometric function (see 13. ; 



F(l/2, -1/2, 1; V 2 ) = -£(7?) 

7T 



(39) 



Similarly, 



dF. 



2) 



ds 



s=0 



F(l/2, 3/2, 2; ^) = ^ (-4) 



k=0 

d 

W) 



v ^ (v 2 ) k (§)*(§)* 



fc=0 



fc! 



(2)a 



F(-l/2, 1/2, 1; ^) = -p- ^E(,,) = - A [E(r?) - K( V )} 



7T rf(?7 2 ) 



7T77 



(40) 



and here we used the formula for the derivative of the hypergeometric function (see |15| , 2.8, (20)), formula (39) again, 
and a formula for the derivative of E(r?) (see JliJ], 13.7, (12)); K(?y) is the complete elliptic integral of the first kind. 
Equations (|38|), (Eoh now provide 



dFs_ 
ds 



_ dF^ 


ofP 


2 


a=0 9s 




s=0 Tf? 



1 + r? 
1 — r/ z 



which coincides with the exact expression in (pL4|) ; the asymptotic formula there for small 6 = 1 — T) is obtained by 
using the expansions of elliptic integrals in the series in the conjugate modulus (see 0, 773.3, 774.3). 
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